home *** CD-ROM | disk | FTP | other *** search
/ Mac Format 1995 June / MacFormat 25.iso / Shareware City / Developers / fortran-to-c-translator-11 / Mac F2C 1.1 / Mac F2C Libraries / libF77 Sources / erf_fctns.c < prev    next >
Text File  |  1995-01-28  |  4KB  |  203 lines
  1. /*
  2.     Error functions required to complete f2c libF77 support library
  3.     
  4.     IMT 4 Jan 94
  5.     
  6.     Modified to match changes in libF77; IMT 30 Nov 94
  7. */
  8.  
  9. #include <math.h>
  10. #include "f2c.h"
  11.  
  12. /* Prototypes for erf functions */
  13.  
  14. double erf( double x );
  15. double erfc( double x );
  16.  
  17.  
  18. /* Prototypes for local functions */
  19.  
  20. static void gcf( double *gammcf, double a, double x, double *gln );
  21. static void gser( double *gamser, double a, double x, double *gln );
  22. static double gammq( double a, double x );
  23. static double gammp( double a, double x );
  24. static double gammln( double xx );
  25.  
  26.  
  27.  
  28. double erf( double x )
  29. {
  30.     return  (x < 0.0) ? -gammp( 0.5, x*x ) : gammp( 0.5, x*x );
  31. }
  32.  
  33. double erfc( double x )
  34. {
  35.     return  (x < 0.0) ? 1.0 + gammp( 0.5, x*x ) : gammq( 0.5, x*x );
  36. }
  37.  
  38.  
  39. /*
  40.     Incomplete Gamma Function  P(a,x)
  41. */
  42.  
  43. static double gammp( double a, double x )
  44. {
  45.     double    gamser, gammcf, gln;
  46.  
  47.     if ( x < 0.0 || a <= 0.0 )
  48.         /* This case can't happen when called by derf() or derfc() */ ;
  49.         
  50.     if (x < (a+1.0))
  51.     {
  52.         gser( &gamser, a, x, &gln );
  53.         return  gamser;
  54.     }
  55.     else
  56.         {
  57.         gcf( &gammcf, a, x, &gln );
  58.         return  1.0 - gammcf;
  59.         }
  60. }
  61.  
  62.  
  63. /*
  64.     Incomplete Gamma Function  Q(a,x)  =  1 - P(a,x)
  65. */
  66.  
  67. static double gammq( double a, double x )
  68. {
  69.     double     gamser, gammcf, gln;
  70.  
  71.     if ( x < 0.0 || a <= 0.0 )
  72.         /* This case can't happen when called by erf() or erfc() */ ;
  73.         
  74.     if ( x < (a + 1.0) )
  75.     {
  76.         gser( &gamser, a, x, &gln );
  77.         return  1.0 - gamser;
  78.     }
  79.     else
  80.     {
  81.         gcf( &gammcf, a, x, &gln );
  82.         return  gammcf;
  83.     }
  84. }
  85.  
  86.  
  87. /*
  88.     Returns following:
  89.         gamser    -    Incomplete gamma function P(a,x) evaluated via series
  90.         gln        -    ln(gamma function(a))
  91. */
  92.  
  93. #define  ITMAX 100
  94. #define  EPS 3.0e-7
  95.  
  96. static void gser( double *gamser, double a, double x, double *gln )
  97. {
  98.     int        n;
  99.     double    sum, del, ap;
  100.  
  101.     *gln = gammln(a);
  102.     if (x <= 0.0)
  103.     {
  104.         if ( x < 0.0 )
  105.             /* ERROR: This case can't happen when called by erf() or erfc() */ ;
  106.         *gamser = 0.0;
  107.         return;
  108.     }
  109.     else
  110.     {
  111.         ap = a;
  112.         del = sum = 1.0/a;
  113.         for ( n = 1; n <= ITMAX; n++ ) 
  114.         {
  115.             ap += 1.0;
  116.             del *= x/ap;
  117.             sum += del;
  118.             if ( fabs(del) < fabs(sum)*EPS )
  119.             {
  120.                 *gamser = sum*exp( -x + a*log(x) - *gln );
  121.                 return;
  122.             }
  123.         }
  124.         /* ERROR: a too large, ITMAX too small.
  125.            This case can't happen when called by derf() or derfc() */ 
  126.         return;
  127.     }
  128. }
  129.  
  130.  
  131. /*
  132.     Returns following:
  133.         gamcf    -    Incomplete gamma function Q(a,x) evaluated via
  134.                     continued fraction
  135.         gln        -    ln(gamma function(a))
  136. */
  137.  
  138. #define  ITMAX 100
  139. #define  EPS 3.0e-7
  140.  
  141. static void gcf( double *gammcf, double a, double x, double *gln )
  142. {
  143.     int        n;
  144.     double    gold = 0.0, g, fac = 1.0, b1 = 1.0,
  145.             b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
  146.  
  147.     *gln = gammln( a );
  148.     a1 = x;
  149.     for ( n = 1; n <= ITMAX; n++ )
  150.     {
  151.         an = (double) n;
  152.         ana = an - a;
  153.         a0 = (a1 + a0*ana)*fac;
  154.         b0 = (b1 + b0*ana)*fac;
  155.         anf = an*fac;
  156.         a1 = x*a0 + anf*a1;
  157.         b1 = x*b0 + anf*b1;
  158.         if ( a1 )
  159.         {
  160.             fac = 1.0/a1;
  161.             g = b1*fac;
  162.             if ( fabs( (g-gold)/g ) < EPS )
  163.             {
  164.                 *gammcf = exp( -x + a*log(x) - *gln )*g;
  165.                 return;
  166.             }
  167.             gold = g;
  168.         }
  169.     }
  170.     /* ERROR: a too large, ITMAX too small.
  171.        This case can't happen when called by derf() or derfc() */ 
  172.     return;
  173. }
  174.  
  175. /*
  176.     Natural log of gamma function.  Full accuracy obtained for xx > 1.
  177.     For 0 < xx < 1, use reflection formula first.
  178. */
  179.  
  180. static double gammln( double xx )
  181. {
  182.     double        x, tmp, ser;
  183.     int            j;
  184.  
  185.     static double    cof[6] = { 76.18009173, -86.505322033, 24.01409822,
  186.                                -1.231739516, 0.120858003e-2, -0.536382e-5 };
  187.  
  188.     x = xx - 1.0;
  189.     tmp = x + 5.5;
  190.     tmp -= (x+0.5)*log(tmp);
  191.     ser = 1.0;
  192.     for (j = 0; j <= 5; j++)
  193.         {
  194.         x += 1.0;
  195.         ser += cof[j]/x;
  196.         }
  197.     return  -tmp+log(2.50662827465*ser);
  198. }
  199.  
  200.  
  201.  
  202.  
  203.